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We employ three numerical methods to explore the motion of low Reynolds number swimmers, 
modeling the hydrodynamic interactions by means of the Oseen tensor approximation, lattice Boltz- 
mann simulations and multiparticle collision dynamics. By applying the methods to a three bead 
linear swimmer, for which exact results are known, we are able to compare and assess the effective- 
ness of the different approaches. We then propose a new class of low Reynolds number swimmers, 
generalized three bead swimmers that can change both the length of their arms and the angle be- 
tween them. Hence we suggest a design for a microstructure capable of moving in three dimensions. 
We discuss multiple bead, linear microstructures and show that they are highly efficient swimmers. 
We then turn to consider the swimming motion of elastic filaments. Using multiparticle collision 
dynamics we show that a driven filament behaves in a qualitatively similar way to the micron-scale 
swimming device recently demonstrated by Dreyfus et al. 

PACS numbers: 



I. INTRODUCTION 

Microscopic and mesoscopic organisms such as bacte- 
ria operate at length scales where swimming motion takes 
place at very low Reynolds number In his 'scallop' 
theorem of microscopic swimming Purcell argued that 
swimming strategies can only be successful in this regime 
if they involve a cyclic and non time reversible motion 
0, 0] ■ The driving of helically shaped bacterial flagella 
by a reversible rotary engine and the beating motion of 
elastic rod-like flagella utilized by eukaryotic cells are ex- 
amples of biological mechanisms which break time re- 
versible invariance, thus allowing microscopic organisms 
to move in a controlled fashion. In an exciting recent de- 
velopment Dreyfus et al. [l[ have demonstrated for the 
first time the controlled swimming motion of a fabricated, 
micrometer size device. 

Several authors have described models of swimmers at 
low Reynolds number. Simple models which comprise 
linked spheres or connected rods that move by chang- 
ing the distances or directions between the components 
are considered in d, H 0, S @1- For one dimensional 
motion analytic results for the swimming velocity and 
efficiency can be obtained. Felderhof @ used the Oseen 
tensor formalism to model microscopic swimmers using 
a bead spring model. Gauger and Stark [l(| used a simi- 
lar method to model the experimental elastic filament of 
Dreyfus et al. Both of these approaches are distinct 
from the one we use, in that the actuation of the beads 
is described in terms of forces, whereas we define swim- 
ming in terms of a predefined shape change. Propulsion 
by a non time reversible pattern of surface distortions 
is addressed in [ill, E • Another possible swimming 
mechanism, mediated by an asymmetric distribution of 
reaction products, is proposed in [Til ]. 

Increasingly, quantitative experiments on bacterial dy- 



namics are appearing in the literature. Transient collec- 
tive motion has been observed in collections of swim- 
ming cells [HI] . Bacteria near solid boundaries have been 
shown to swim in circles [1 6j j and those near an obstacle 
to reverse their swimming direction [17|. Experiments 
have been performed to determine the dependence of the 
chemotactic response of Dictyostelium discoideum cells 
swimming in a concentration gradient [181 ]. 

Although simple, analytical models can provide con- 
siderable help in understanding these results there are 
many new features inherent or accessible in real bio- 
logical systems that remain to be explored. These in- 
clude more complicated swimming mechanisms, interac- 
tions between densely packed swimmers and the effect of 
boundaries and obstacles. There will increasingly be a 
need to develop numerical methods to probe the behav- 
ior of more complicated structures and situations where 
analytic approaches become intractable. 

In numerical approaches published so far Hernandez- 
Ortiz et al. have considered the swimming motion of a 
collection of force dipoles 0. They observe the large 
scale coherent vortex motion that has been seen in ex- 
periments. Ramachandran et al. have described swim- 
mers, modeled as force dipoles, interacting with a fluid 
described by a lattice Boltzmann algorithm [20(. Work 
on swimmers propelled by a flexible filament modeled hy- 
drodynamics through an anisotropic friction coefficient 

HHI2I- 

Our first aim in this paper is to explore new ways in 
which simulation methods can be applied to the motion of 
swimmers in a low Reynolds number solvent. We model 
the hydrodynamic interactions by using the Oseen ten- 
sor approximation [231 ]. a lattice Boltzmann algorithm 
[2~il l25j and multiparticle collision dynamics 26]. The 
approaches are validated by solving the equations of mo- 
tion for a linear three bead swimmer where an analytic 
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solution is available for comparison Q. We discuss the 
relative merits and demerits of the three approaches. 

Secondly, we extend the linear model to more general 
three bead microstructures. We use an Oseen tensor ap- 
proach to demonstrate that they can move in a controlled 
fashion in three dimensions by changing both the length 
of and the angle between their arms and we discuss the 
efficiencies of the various swimming modes. We also show 
that multibead linear structures are highly efficient swim- 
mers. 

We then consider the swimming motion of driven elas- 
tic filaments. Our model, solved using multiparticle col- 
lision dynamics, mirrors the behavior of the swimming 
device introduced by Dreyfus et al. []J. In the conclu- 
sion, we consider future directions in which the modeling 
approaches might be useful in understanding the swim- 
ming of bacteria and of fabricated microstructures. 



II. NUMERICAL APPROACHES 

The swimmers we consider are composed of N spheres, 
of fixed radius R and with positions given by rj, where 
i = 1...N. The spheres are linked by rods that are suffi- 
ciently thin to neglect any hydrodynamic effect. Internal 
forces and torques act to change the lengths and/or an- 
gles between the rods, causing the swimmer to change 
shape. These shape changes, when coupled to the fluid, 
lead to directed motion. We now describe three different 
numerical approaches used to simulate the fluid. 



A. Oseen tensor 

The Oseen tensor allows us to consider the hydrody- 
namic interaction, in the limit of zero Reynolds number, 
between spheres that are spaced far apart (i.e. at dis- 
tances significantly larger than their radii, R) [23|]. A 
sphere pushed by a force will move, and so set up a flow 
field. Any surrounding spheres will be advected with 
the resulting local fluid velocity. Furthermore, since the 
Reynolds number is very low, the time taken to set up the 
flow fields is much smaller than that needed for a given 
sphere to move a significant fraction of R [27[. There- 
fore, the hydrodynamic interaction can be approximated 
as instantaneous. Since the Stokes equation (the Navier- 
Stokes equation without the inertial term, as is appro- 
priate in the low Reynolds number limit) is linear, the 
velocity fields produced by each of the spheres simply 
add up. This allows us to write 

N 




FIG. 1: (a) The swimmer at time t. (b) The swimmer shape 
defined at time t + St. The lengths of the two links connecting 
sphere 3 have decreased, (c) The shape at t + St is translated 
by a vector Ar and rotated by an angle A9 around its center 
of mass. These operations preserve the link lengths, thus 
leaving the shape unchanged. The parameters Ar and A8 
are chosen to improve upon the accuracy of the constraints 
in equation Q. This procedure is performed iteratively until 
the necessary accuracy is achieved. 

where v$ is the velocity of sphere i and Pj is the force on 
sphere j. The Oseen tensor is (23| 

lJ 1 5-4— rfl+r^) otherwise, K ' 

K 87r?7|ry| V Nil 2 J ' 

where r\ is the fluid viscosity, I is the identity matrix, 
and Yij = Tj — Yi is the vector between spheres i and j. 
For a swimmer, the forces Fj are not external, but rather 
internal forces mediated through the links that connect 
the spheres. They are subject to the constraints 

N N 

^F, = 0, ^F.xr^O (3) 

i=i i=i 

which state that no external forces or torques act on the 
swimmer. 

The swimming motion is defined as a periodic shape 
change and, from this information, the algorithm must 
determine the trajectory of the swimmer through the 
fluid. To illustrate how this works we begin by con- 
sidering a swimmer whose motion is confined to a two 
dimensional plane. Figure [T] shows the procedure for the 
case N = 3. At a given time t, the position of the spheres 
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Yi are known, as shown in figure QJa). The new shape of 
the swimmer at the next time step t+St is shown in figure 
[ljb). We have chosen, for illustrative purposes, a swim- 
ming step where the lengths of the two links connecting 
sphere 3 have decreased. 

The shape of the swimmer in figure W[J°) is defined 
through the three quantities \r' 12 \, l r 23l; an d l r 3il- How- 
ever, this information does not determine the absolute 
positions of the spheres. To find these, it is necessary 
to enforce the conservation conditions stated in equation 
([3]). This is performed in the following iterative manner. 

We take the first approximation for Ti{t + St) to be rj. 
This does not, in general, obey equation fl3|), as will be 
apparent below. Our aim is to move the spheres in such a 
way as to successively improve the accuracy of equation 
([3]). To do this, we consider translating the swimmer by a 
vector Ar = (Ax,Ay) T , and rotating it about its center 
of mass by an angle A9, as illustrated in Figure HJc). 
Note that these operations do not change the shape of 
the swimmer. In doing this, we introduce new position 
coordinates r" defined by 



used in equation ([4]) to obtain a better approximation, 
r", to the sphere positions at time t + St. Repeating the 
procedure gives rapid convergence to the correct solu- 
tion. Typically a single iteration improves the accuracy 
of equation ([3]) by a factor of ~ 10. 

The approach is easily generalized to a swimmer mov- 
ing in three dimensions. In this case there are three dis- 
placement parameters Ax, Ay, Az, and three rotation 
parameters A9 xy A9 yz , A9 ZX (where, for example, A9 xy 
is a rotation in the x-y plane). These six unknowns can 
be determined using the six constraints in equation ([3]). 
The method is the same as above, except equations 
[5]) now contain the corresponding extra terms, and A in 
equation {3 becomes a six by six matrix. 

The computation for one time step scales as TV 3 . For 
low to moderate values of N the method is very fast. 
However, for N > 100, it quickly becomes very compu- 
tationally intensive. 



Lattice-Boltzmann 



r + Ar - 

A cm ' " A 



(4) 



where R is a clockwise rotation matrix around an angle 
A9 and r' cm = Y^Li V 'i/N is the center of mass position 
of the swimmer. The unknown displacement Ar and 
angle A9 are chosen to improve upon the accuracy of 
the constraints in equation ([3]). To calculate them we 
note that the velocity of sphere i can be related to its 
displacement over time St by 



St 
r' - r, 



Ar + A0K-rUr? 



St 



(5) 



The unit vector if lies in the direction a given point 
moves under an infinitesimal rotation (which can be cal- 
culated by rotating the vector — r' cm by 90° clockwise 
and normalizing it). 

Substituting expression © into equation ((T|), and nu- 
merically inverting the resulting matrix equation using 
Gaussian elimination, we obtain an expression for the 
forces acting on each sphere of the form 



a.i + b 4 Az + c t Ay + d t A9 



(6) 



where a^, b,, Cj, and are constant vectors. Using 
this, the three constraints in equation ([3]) can, after some 
rearranging, be written 



(7) 



where the matrix A and the column vector e are con- 
stants. Finally, this matrix equation can be inverted to 
obtain the values of Ax, Ay, and AO which are then 





The lattice Boltzmann algorithm is now a widely used 
mesoscopic modeling technique for simulating the behav- 
ior of complex fluids [24|, |25| ■ The method consists of an 
evolution equation for a mass density distribution func- 
tion f k (s,t), which can be considered as a simplified, dis- 
cretized version of Boltzmann's transport equation. The 
distribution function is defined at positions, s, which lie 
on a cubic lattice with a distance Ss between nearest 
neighbor points. Its value is updated simultaneously and 
discretely in time, with time step St. We define c = Ss/ St. 
The subscript k denotes a particular velocity direction . 
The velocity vectors must be chosen such that e k St lies 
between lattice sites. In this study all simulations are 
performed in three dimension using a 15 velocity model. 
This has a zero velocity vector eo = (0, 0, 0), six nearest 
neighbor velocity vectors ei_ 6 in the directions (±c, 0, 0), 
(0, ±c, 0), and (0, 0, ±c), and eight velocity vectors e 7 _i 4 
in the diagonal directions (±c, ±c, ±c). From //. we can 
calculate the mass density p and momentum density pu: 



pu a 



(8) 



Evolution in time is given by 



A(s + e k St, t + St) = f k (s, t) - - f k (s) - f^(s) (9) 

T L J 

where we use the Bhatnagar-Gross-Krook approxima- 
tion, which uses a single parameter r to determine the 
rate of relaxation toward equilibrium. A suitable choice 
for the equilibrium distribution function is 



fl q = pw k [l + 



3e ka u a 9(e ka u a ) 2 3-u^ 



2c 4 



2c 2 



where the weight factors are wq = g, wis = jj) 



(10) 



and 
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W7-14 = Note that this distribution satisfies 

Yl -fk q = P> fftkci = pu a , (11) 

k k 

such that mass and momentum are conserved in time. 
This can be seen by summing the zeroth and first velocity 
moments of equation ([9]), and using the relations in ([8]), 
Applying the Chapman-Enskog expansion to the lat- 
tice Boltzmann equation ([9]) [24( gives the continuity 
equation for the total density, 

d tP + d Q (pu a ) =0 (12) 

and the Navier-Stokes equation for the fluid momentum, 

d t {pu a ) + dpipuaup) = -d a (pc 2 s ) + d/s {iypdpu a ) (13) 

where the kinematic viscosity is 

and the speed of sound is c s — cj \/3. At low Reynolds 
numbers the fluid is, to a very good approximation, in- 
compressible, i.e. V.u = 0. (In the simulations, the 
difference between the maximum and minimum density 
within the system was found to be no more that 0.04% 
of the average density.) 

For simplicity, we do not explicitly consider a solid- 
fluid interface but define the swimmer to comprise of the 
fluid reg ion within the spheres which make up the swim- 
mer [28l | . The swimmer- fluid interaction is incorporated 
into the lattice Boltzmann algorithm after finding the 
fluid velocities in equation (8), but before calculating the 
equilibrium distributions in equation (10). This interac- 
tion is generated in three stages. Firstly, the total linear 
and angular momentum of the swimmer are calculated: 

P = Pi u i ' L = S J x Pi u i ' ( 15 ) 

3 3 

where the sum j runs over all the lattice sites contained 
within the swimmer. Secondly, the new positions of the 
spheres are calculated. This procedure is analogous to 
that described for the Oseen tensor method in Section 
III Al In this, we know the positions of the spheres at time 
t and the swimmer shape at time t + St. The algorithm 
works out how this new shape is oriented with respect to 
the original, such that linear and angular momentum are 
conserved, i.e. 

^ rriiVi = P, Y ViX nilVi = L ' ( 16 ) 

i i 

where mi — • pj and 

_ r i (t + St)-r i (t) 

St (U) 

are the mass and velocity of sphere i, respectively. 
Thirdly, the motion of the swimmer is coupled back to the 



fluid. Lattice sites within a given sphere are set to the ve- 
locity of that sphere, i.e. Uj = v^, Sj 6 Sphere i. These 
updated velocities are then used in calculating the equi- 
librium distributions (|10p . Through repeated iteration 
of the lattice Boltzmann equation (9), the fluid within 
and immediately adjacent to the spheres is strongly cou- 
pled to move with the same velocity as the spheres, thus 
giving the correct boundary condition. 

To avoid unwanted lattice effects, the radius R of each 
sphere must be sufficiently large to accurately resolve its 
shape on the cubic grid. In this study we choose R = 3Ss, 
such that each sphere contains approximately 113 lattice 
sites. Note that in this procedure we neglect the effect of 
rotation on the spheres, assuming that all parts of a given 
sphere travel at the same velocity. This assumption is 
justified because the hydrodynamic interactions between 
rotating spheres is rather weak. (The Oseen tensor ([2]) 
decays as r _1 , whereas the flow field around a rotating 
sphere decays at a much faster rate of r~ 3 .) In the case 
of modeling more than one swimmer, it is necessary to 
add hard core potentials between spheres to prevent them 
from overlapping. 

C. Multiparticle Collision Dynamics 

An alternative mesoscale approach, which solves the 
equations of fluctuating hydrodynamics is the multipar- 
ticle collision dynamics (also known as stochastic rota- 
tion dynamics) algorithm introduced by Malevanets and 
Kapral (26|. The fluid is represented by a large number 
of point-like particles of mass m, with position r/j(i), and 
velocity Vk(t) at time t, where k is the particle index. The 
particles move in continuous space, and at discrete time 
intervals, St. Particle positions are updated according to 

r k (t + St)=r k (t)+v k (t)5t. (18) 

At each time step the particles also undergo a multiparti- 
cle collision that locally conserves mass, momentum, and 
energy. To perform the collision, the simulation box is 
divided into a grid of cubic cells, with sides of length a. 
The average number of particles per cell will be denoted 
by 7. The velocities of the particles in each cell are then 
rotated about the center of mass velocity of the cell, v cm 

v fc (t + St) = v cm (i) + R [vfe(i) - Vcm (t)] . (19) 

R is a rotation matrix through a fixed angle, a, about 
an axis that is generated randomly for each cell in the 
simulation at each time step. The position of the cubic 
grid is chosen randomly at each time step - this leads 
to substantially improved Galilean invariance in the al- 
gorithm [29l ]. In the continuum limit, the multiparticle 
collision dynamics algorithm recovers the thermohydro- 
dynamic equations of motion and thus acts as a Navier- 
Stokes solver. Conveniently, the dependence of the trans- 
port coefficients of the fluid on the simulation parameters 
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is known analytically [3(J, [H|. Thus, it is a relatively sim- 
ple task to choose values that result in a low Reynolds 
number fluid. 

We couple a swimmer to the multiparticle collision dy- 
namics solvent by considering it to be composed of a 
number of particles, and including these solute particles 
in the solvent collision step (fH))) . In this way the swim- 
ming microstructure can exchange momentum with the 
solvent. In general, the equations of motion of the mi- 
crostructure are solved using a velocity Verlet molecular 
dynamics algorithm. Precise details of how specific swim- 
mers are dealt with using this approach are given in later 
sections of the paper. 



III. LINEAR THREE SPHERE SWIMMER 

Recently, Najafi and Golestanian pj proposed a 
one-dimensional swimmer comprising three connected 
spheres. Their model is perhaps the simplest example 
of a controlled, cyclic motion that breaks time reversibil- 
ity. The swimmer consists of a central sphere that is 
connected to two other spheres by arms that are sep- 
arated by an angle of 180°, are of negligible thickness, 
and whose length can be changed by, for example, the 
action of engines located on the central sphere. The 
microstructure moves by shortening and extending the 
lengths of the arms in a periodic and time irreversible 
manner, as shown in Figure O The relevant parameters 
for this model are D, the distance between the central 
sphere and an outer sphere at the maximum arm length, 
e, the distance the arm shortens, W, the speed at which 
the arms change their lengths, and R, the radius of the 
spheres. The result of this cyclic, time irreversible mo- 
tion is a net translation of the swimmer along the line 
linking the three spheres; we define A as the distance 
the swimmer translates in one complete cycle. 




x 

► 



FIG. 2: The four step, cyclic motion of the linear three sphere 
swimmer Q. 



undergoes small deformations, is 



A. Analytic Theory 

Because of the simplicity of the shape deformations of 
the linear three sphere swimmer, it is possible to calculate 
analytically the total net displacement, A, of the swim- 
mer during each complete cycle of its motion, in the limit 
of e <C D and R <C D. We summarize the argument of 
Najafi and Golestanian [tJ and correct their formula for 
the displacement of the swimmer, which we shall need 
for comparison to the numerical results. 

Each of the four steps of the stroke can, by a simple 
transformation, be converted into one particular auxil- 
iary stroke. In the auxiliary stroke one arm has a fixed 
length, S, where 6 is either D or D — e, and the other arm 
changes length from D to D — e. We choose the a;- axis 
to be parallel to the line linking the spheres and directed 
away from sphere 2 (see Figure [5]) • During the auxiliary 
stroke v\ — v% and W = V2 — V\, and thus the veloc- 
ity of the middle sphere, in the limit that the swimmer 



vi(5) 



-W (H u - H 23 - H 12 ) 



(3-ffn 
_W 



2 (H12 + Hi 3 
R 

~ 2{D - Wt) 



H23)) 
R 

~S ~ 2(5" 



R 



D 



Wt) 
(20) 

ignoring terms of order (R/D) 2 and greater. The ele- 
ments of the Oseen tensor for each pair of spheres follow 
from equation ([2]). Integrating (|20|) gives the displace- 
ment over the auxiliary stroke, 
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(21) 



This can then be used to calculate the total displacement 
after the four step cycle, to second order in e/D, as 
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We note here that this expression differs from that 
given by Najafi and Golestanian [3] who reported 



A = 2.8R ( — 



(23) 



In equation (|23[) the displacement, A, is proportional to 
(e/D) 3 . If one considers the transformation e — > x = — e 
and -D — > G = D — e, this corresponds to a swimmer 
undergoing exactly the same continuous motion as that 
shown in Figure [H only with the swimming stroke begin- 
ning at the third step in the cycle. Thus, the swimmer 
must move in the same direction. However, equation l|23j) 
suggests that the swimming direction is reversed under 
this transformation, which is clearly incorrect. 

In the following section, we summarize the results of 
numerical simulation studies of the motion of the linear 
three sphere swimmer at low Reynolds number in order 
to validate and compare the use of these approaches in 
the study of swimming microstructures. 



B. Numerical Results 

Figure[3jgives the results for a single linear three sphere 
swimmer, using each of the three methods presented in 
Section HU This graph shows how the total displacement 
of the swimmer over one swimming cycle, A, varies as a 
function of the amplitude of the stroke, e. The parame- 
ters used were D = 25 and R = 3 for the Oseen tensor 
and lattice Boltzmann approaches and D = 3.0a for the 
multiparticle collision dynamics simulations. 

The solid line in Figure[3jis obtained by directly solving 
the Oseen tensor interaction between the spheres, as out- 
lined in Section Hi Al The dashed line shows the theoret- 
ical curve, equation (|22[) . which is correct to third order 
in e/D. These two curves converge in the limit of small 
e/D, as expected. The dot-dashed line is equation (|2"3l . 
the expression proposed by Najafi and Golestanian 
This appears to give good agreement for larger values of 
e/D. This is misleading, however, as in the limit of small 
e/D it does not converge to the theoretical solution (this 
is seen most clearly within the inset), and it should not 
be valid at higher values of e/D due to the assumptions 
made in the derivation. 

The lattice Boltzmann simulations were performed us- 
ing a lattice of size L x = 200, L y — 100, and L z = 100, 
with periodic boundary conditions. Initially, the swim- 
mer was placed in the middle of the box, aligned parallel 
to the x direction. The relaxation parameter was chosen 
to be t = 1. The simulations were run for one complete 
swimming cycle, which corresponded to t max = 102400<5t 
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FIG. 3: The shift per cycle of the linear three sphere swim- 
mer, A, as a function of the sphere displacement amplitude, 
e. The inset shows a magnified view of the region of the graph 
below e/D — 0.2. The parameters D = 25 and R — 3 were 
used in the Oseen tensor and lattice Boltzmann approaches. 
In the multiparticle collision dynamics simulations we used 
D = 3.0a. The solid line was obtained by numerically solv- 
ing the Oseen tensor equation, outlined in section Hi Al The 
crosses mark results obtained from lattice Boltzmann simu- 
lations, presented in section III Bl The error bars show the 
distribution of results using multiparticle collision dynamics 
from section III CI The dashed line is the theoretical expres- 
sion given in equation (|22[) and the dot-dashed line is the ex- 
pression proposed by Najafi et al. Q reproduced in equation 

us. 



time steps. The maximum speed of spheres is approx- 
imated by ^e/t max . Using this together with R, which 
gives a characteristic length scale for the problem, the 
Reynolds number can be expressed as Re = 4eR/iyt max . 
The largest displacement used (e = 19) gives Re — 0.013. 
This was checked to be sufficiently low by running a lim- 
ited number of simulations using t max = 2048005i time 
steps, and finding that these results agreed to within 1%. 
Furthermore, we checked that finite size effects were not 
important by performing simulations using a lattice of 
size L x = 300, L y = 150, and L z — 150, with results 
again agreeing within this tolerance. 

The results from the lattice Boltzmann simulations are 
denoted by the crosses in Figure [31 They agree well with 
the full Oseen tensor result (solid curve) at small e/D and 
deviate as e/D gets increasingly large. This is because 
the Oseen tensor approximation that the spheres are far 
apart breaks down in this limit (the spheres intersect each 
other if e/D > 0.76 for the parameters used). 

Simulations using multiparticle collision dynamics are 
intrinsically noisy. If we simply use a molecular dynamics 
algorithm to solve the equations of motion for the swim- 



mer and include the particles that make up the swimmer 
in the solvent collision step then, without further correc- 
tion, the microstructure will rotate, and not stay aligned 
along one particular axis. Although this behavior would 
be realistic for a nanoscalc microstructure in a solvent, 
it does not easily allow for an accurate comparison of 
swimming speed with theory. To constrain the motion 
to one dimension, the transverse velocities of the three 
particles that comprise the swimmer were adjusted to the 
average velocity of the particles after each collision step. 
Changes in the arm lengths were undertaken by adding 
an extra velocity to each particle, such that the total 
momentum of the microstructure remained unchanged 
during the arm length change. 

For the multiparticle collision dynamics solvent we 
used the following parameters: Particle temperature 
kT = 0.005, time step St = 0.01, cell size a = 1.0, 
rotation angle a — 135°, average number of particles 
per cell 7 = 10, and particle mass m — 10. These 
result in a Reynolds number for the microstructure of 
~ 10~ 5 . These parameters both ensure a low Reynolds 
number and minimize fluctuations in the solvent with a 
high Schmidt number. In our simulations, we employed 
a simulation box of dimensions 30a x 8a x 8a with peri- 
odic boundary conditions and checked for finite size ef- 
fects. For the swimmer we used a mass of 5m for each 
sphere. Due to the nature of the swimmer-solvent in- 
teraction in our implementation of the multiparticle col- 
lision dynamics algorithm, it is difficult to define the 
effective hydrodynamic radii of the spheres. However, 
comparison with the Oseen tensor and lattice Boltzmann 
simulation results suggests the parameters used lead to 
R/D ~ 0.12. The simulations were conducted for a total 
time of t m ax = 2.72 x lO 5 ^ time steps, and one period of 
motion took 6.8 x 10 3 <5t. The period must be sufficiently 
long to allow the solvent to couple with the swimmer. 20 
runs were performed for each parameter set and the re- 
sults are denoted by the error bars in Figure G2 which are 
spread one standard error on either side of the average of 
the 20 runs. The results are compatible with, but much 
less precise than, those obtained by the methods without 
intrinsic fluctuations. 



By changing the parameters, D and e, it is not only 
possible to change the swimmer displacement, A, but 
also the efficiency of the swimmer. We define this effi- 
ciency to be the energy required for an external force to 
move the individual spheres by the distance A in a time 
P divided by the work done by the swimmer in perform- 
ing the corresponding cyclic shape change: 



Efficiency : 



6iri]NA 2 /P 
Eti Jo 'Fi.Vidt 



(24) 




Using the Oseen tensor approach, the forces on the 
spheres, Fj, are calculated through equation ©, so this 
quantity is easily obtainable. The line in figure HJa) 



FIG. 4: The percentage efficiency of a linear three sphere 
swimmer as a function of (a) the sphere displacement ampli- 
tude, e, for fixed D = 5R, and (b) the maximum arm length, 
D, for fixed minimum separation D — e = 4_R. 



shows how the swimming efficiency is changed as a func- 
tion of e, whilst keeping D = 5R fixed. Note that this 
curve terminates when e = D — 2R, since beyond this 
the spheres unphysically overlap at some point within the 
swimming cycle. Furthermore, the Oseen tensor method 
tends to overestimate the efficiency when the spheres get 
close, because it does not include the viscous dissipation 
resulting from lubrication effects. Figure Ufa) illustrates 
a general feature of swimmers, namely that small ampli- 
tude motion results in inefficient swimming. The curve 
in figure HJb) shows the efficiency against -D, assuming a 
fixed mimumum sphere separation of D — e = AR. We 
find that the swimmer becomes increasingly efficient as 
D increases, approaching a limit of around 8%. 

We now summarize the relative advantages and disad- 
vantages of the three numerical methods. The Oseen ten- 
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sor approach is particularly advantageous for swimmers 
comprising small numbers of spheres because it is com- 
putationally very fast (simulations taking only minutes 
instead of days using lattice Boltzmann or multiparticle 
collision dynamics). This is primarily because it is not 
necessary to explicitly solve a set of fluid dynamics equa- 
tions. However the simulation time scales as iV 3 so lattice 
Boltzmann becomes more efficient for large numbers of 
swimmers. 

The Oseen tensor formalism is also limited because it 
assumes that the interacting spheres are spaced far apart, 
compared to their radii. This means, for instance, that it 
would not be appropriate to use this method to study the 
movement of a swimmer close to a wall or in a confined 
geometry. On the other hand, the lattice Boltzmann al- 
gorithm can address these problems, providing an exact 
solution to any fluid flow problem given sufficient resolu- 
tion. In practice it is limited by computational power. To 
avoid spurious lattice effects, the sphere radius must be 
significantly larger than the lattice size, necessitating the 
need for large systems. Moreover, to obtain a sufficiently 
low Reynolds number, the cyclic swimming motion must 
be performed over a great number of time steps, further 
increasing the computational burden. 

Multiparticle collision dynamics is advantageous be- 
cause it is unconditionally stable, unlike the lattice Boltz- 
mann method, and this allows somewhat lower Reynolds 
numbers to be obtained more easily. One can also use a 
molecular dynamics approach to treat the microstructure 
allowing for great flexibility in the swimmer models it is 
possible to consider. 

Additionally, as this method inherently contains noise 
it will be appropriate for studying very small scale struc- 
tures, for which Brownian fluctuations are important. 
However, if the fluctuations are unphysical the noise is 
undesirable, necessitating long time averages. 

For the remainder of this paper, we concentrate on us- 
ing the Oseen tensor and multiparticle collision dynamics 
approaches to study more complex swimming motions. 



IV. GENERALIZED THREE SPHERE 
SWIMMERS 

The swimmer described in Section IIIII is constrained 
to move in one dimension along its axis. Using the same 
basic elements one can design a number of other three 
sphere swimmers that can move their individual compo- 
nents and centers of mass in two or three dimensions. To 
extend the original design we allow the angle between the 
two arms of the swimmer to change @ . When the change 
in angle takes place, the spheres may move radially with 
the arm lengths constant, or the arm lengths may change 
at the same time, thus allowing for a number of differ- 
ent motions. In Figure [5] we show one of many possible 
alternative schemes of motion for swimmers of this type, 
which we will refer to as generalized three sphere swim- 
mers. 




X 



FIG. 5: Alternative cyclic motion for a three sphere swimmer, 
allowing the microstructure to translate in two dimensions. 
Possible extensions to this scheme include allowing the angle 
a to change as the arm lengths change. 
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FIG. 6: The position shift over one cycle, A, as a function of 
swimmer angle a for the swimmer defined in Figure [5] Other 
parameters were D = 57? and e = 2R. 

We first concentrate on the motion shown in Figure 
and the case where the arms rotate at fixed length. If the 
swimmer employs a symmetric motion with arm lengths 
D\2 = D13 and s\2 — £13, the net translation of the 
swimmer is along the x-axis (defined in the figure). Thus, 
the microstructure remains a one dimensional swimmer 
while its individual elements are moving in two dimen- 
sions. Interestingly, the direction that the microstructure 
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FIG. 7: The angular change in the orientation of the gen- 
eralized three sphere swimmer over one complete swimming 
cycle, AO, as a function of the swimmer displacement ampli- 
tude, £13. Other parameters were fixed to be £12 = 3R and 
D12 — D\3 = 6R. The inset shows the movement of the cen- 
ter of mass of the swimmer at a fixed point in the cycle, over 
the course of 9 swimming cycles. 



translates varies with a, for fixed D and e, as shown in 
Figure O For D = 5i? and e = 2R the outer two spheres 
do not intersect for a > 39°, and there is a transition 
from backward motion to forward motion at a = 77°. 
Intuitively it is not obvious what causes this reversal of 
direction. We found that the maximum efficiency of this 
swimmer (with efficiency defined by equation (|24[>) . oc- 
curs at a = 138°. This corresponds to an efficiency of 
1.8%, which is considerably less than the ~ 8% found for 
the linear three sphere swimmer in the previous section. 
It seems plausible that the linear three sphere swimmer 
is the most efficient three sphere swimmer possible. 

If one instead defines an asymmetric motion, with 
D12 7^ D13 and/or £12 7^ £13, the microstructure will 
move its center of mass in two dimensions. The swim- 
mer rotates and translates, as shown for D12 = D13 = 6, 
£12 = 3 and £13 as variable in Figure [71 As the difference 
between £12 and £13 becomes greater, the angular rota- 
tion about the center of mass of the microstructure for 
each swimming cycle increases. Thus, it is a relatively 
easy step to imagine a manufactured device that can 
switch between symmetric and asymmetric cyclic mo- 
tions, perhaps in response to an external stimulus or ex- 
perimental condition, to enable movement in a controlled 
fashion in two dimensions. A three sphere swimmer, that 
can change the angle between its two arms, could sim- 
ply adopt the efficient swimming motion of Section IIII1 
where a = 180°, to move in a straight line, and then vary 
a to adopt a structure allowing it to turn. 

To extend the movement of such a microstructure to 
three dimensions, the angle between the arms could be 



changed along another plane, perpendicular to the orig- 
inal one. One strategy to do this might be to em- 
ploy a double-jointed structure where the angle between 
the arms could be changed from a = 180° in either of 
two perpendicular planes. It is of interest that this mi- 
crostructure is the first low Reynolds number swimmer to 
be proposed theoretically that can move in a controlled 
fashion in three dimensions, without employing numer- 
ous one dimensional swimming devices placed perpendic- 
ular to each other. 



V. EXTENDED, LINEAR, ONE DIMENSIONAL 
SWIMMERS 

To extend the linear three sphere swimmer of Section 
IIIII one can simply add more spheres to the microstruc- 
ture, the simplest extension being the four sphere case. 
For the four sphere microstructure shown in Figure [5J), 
we analyzed through Oseen tensor based numerical sim- 
ulations all possible cyclic motions that are made up of 
a discrete number of steps. At the end of each step the 
distance between neighboring spheres is either D (an ex- 
tended rod) or D — e (a contracted rod). If a rod changes 
length during a step then it does so at a constant speed 
W. In this analysis we allowed for more than one rod 
length changing simultaneously. 

Out of this subset of possible swimmers the swimming 
strategy shown in Figure [8] is the most efficient. This 
optimal swimming strategy proceeds as follows: Starting 
from the fully extended conformation, i), the distance be- 
tween spheres 1 and 2 is reduced to D — e, ii). In the next 
two steps the distance between spheres 2 and 3 is reduced 
to D — £, iii), and the distance between spheres 3 and 4 
is reduced to D — e, iv), the fully contracted conforma- 
tion. The microstructure then sequentially extends, first 
by extending the distance between spheres I and 2 to D, 
v), then extending the distance between spheres 2 and 3 
to D, vi), and finally by extending the distance between 
spheres 3 and 4 to D, completing the cycle and taking 
the conformation back to the original starting configura- 
tion i). For the case D = 5R and e = 2R, the four sphere 
microstructure using the optimal swimming strategy has 
a swimming efficiency of 6.9 % compared to 4.5 % for the 
three sphere swimmer using the same parameters. 

Extending to the five sphere case, we analyzed all 
possible cyclic swimming strategies using Oseen tensor 
based numerical simulations, and found that the anal- 
ogous swimming strategy to that in Figure [8] was opti- 
mal, resulting in a swimming efficiency of 8.8 % for the 
D = 5i? and e = 2R case. We studied this optimal swim- 
ming strategy for microstructures with up to 200 spheres 
and the swimming efficiency as a function of the num- 
ber of spheres is shown in Figure [5J The curve indicates 
a logarithmic growth in the efficiency as a function of 
sphere number. This would imply that for a sufficiently 
large number of spheres the efficiency will go above one. 
Although counter intuitive, from our definition of effi- 
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FIG. 8: The most efficient cyclic swimming strategy for an 
extended, linear, four sphere microstructure. 
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FIG. 9: Percentage efficiency against the number of spheres 
for a microstructure adopting the swimming strategy depicted 
in Figure HI 



ciency in equation (|24[) . this is perfectly possible, and 
does not violate any physical principles. However, as the 
number of spheres is increased, the physical size of the 
spheres and rods would have to be made proportionately 
smaller, to ensure that finite Reynolds number effects do 
not become important. 



VI. FILAMENT SWIMMERS USING 
MULTIPARTICLE COLLISION DYNAMICS 

As an example of using multiparticle collision dynam- 
ics to investigate a more complicated swimmer we con- 
sider the motion of a filament modeled as beads con- 
nected by springs and interacting through Lennard- Jones 
forces. The filament is driven by a sinusoidally oscillat- 
ing force applied at one end. The model was motivated 
by the man-made microscopic swimmer of Dreyfus et al. 
where a red blood cell is attached to a filament consist- 
ing of superparamagnetic colloids, that are connected to 
each other using DNA Q. Two magnetic fields are used 
experimentally, one to align the filament and the other 
to actuate one end of it in a sinusoidal manner. This 
actuation results in a series of waves, originating at the 
tail of the filament, propagating towards the red blood 
cell at the head. Because the perpendicular and parallel 
friction coefficients of the microstructure are not equal 
(Ci/C|| ~ 2), a net translation, in the opposite direction 
to the propagation of the wave, occurs along the align- 
ment direction of the filament. 

In the multiparticle collision dynamics simulation we 
model the filament as a number of Lennard-Jones beads, 
representing the superparamagnetic colloids, connected 
to each other by FENE springs, representing the DNA 
linkers. Instead of a magnetic field, we simply apply an 
equal and opposite force to each end of the filament to 
align it along an axis, and apply a sinusoidal actuating 
force, perpendicular to the aligning forces, to one end of 
the filament (see Figure fTOfa)). 

Lowe [32] proposed a dimensionless parameter to char- 
acterize naturally flexible filaments, the 'sperm number', 
defined as 



L/ 



1/4 



(25) 



where L is the length, k is the bending rigidity, and £j_ 
is the perpendicular friction coefficient for the filament, 
and to is the angular frequency of the actuation or driv- 
ing. For our model filament in a multiparticle collision 
dynamics solvent we calculate the bending rigidity from 
the change in energy of the structure as a function of the 
curvature of the filament ■ We determine the friction 
coefficients by applying a known force to each bead, in 
a direction perpendicular or parallel to the alignment, 
and measuring the resulting velocity of the microstruc- 
ture (with no actuation). We can then measure how the 
velocity of the microstructure depends on S p by changing 
the angular frequency of the actuation. 
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FIG. 10: (a) Model filament used in the multiparticle collision 
dynamics simulations. The filament is composed of Lennard- 
Jones beads that are connected by FENE springs. An aligning 
force, F a iig n , is applied to both ends of the filament, and a si- 
nusoidal force, Factuate = F ma x sin(w£), is used to actuate one 
end of the microstructure, leading to wave propagation from 
right to left in the diagram, (b) Scaled swimming velocity 
as a function of the sperm number for the flexible filament. 
To obtain this plot, F max is varied for different frequencies of 
actuation such that the integral of F ac tuate over half of the 
period is equal. 



For the multiparticle collision dynamics solvent we use 
the following parameters: kT — 1.0, St = 1.0, a = 1.0, 
a = 120°, 7 = 5, and m = 4 resulting in a Reynolds num- 
ber for the microstructure of ~ 10~ 2 . For the microstruc- 
ture we use a mass of 4m for each bead and a molecular 
dynamics time step of 0.002St. The average distance be- 
tween the centers of mass of each bead is « 1.0a, and we 
use 10 beads to represent the filament. The simulations 
are conducted over t max = 200000(5i solvent time steps, 
and we average the results over 10 runs for each data 
point. 

In Figure [TO] we show the swimming velocity of the 
microstructure, scaled by Lui, as a function of S p . As 
predicted by theory for naturally flexible filaments [33j ] . 
we observe a maximum in the scaled swimming velocity 
of the filament between the high (dominated by viscous 
friction) and low (dominated by internal elasticity) S p 
regimes. This reproduces the behavior demonstrated for 
the man-made swimmer of Dreyfus et al. [l[ , and we also 
observe very similar scaled speeds in our simulations to 
those found experimentally. 

To verify that the filament in the simulations swims 
through the mechanism of wave propagation, we also 
modeled a 2 bead microstructure under the influence 
of aligning and actuating forces. We find that this mi- 



crostructure does not swim as it is impossible for a 2 bead 
filament to move in a time irreversible fashion. 



VII. CONCLUSIONS 

In this paper we have described three different meth- 
ods for simulating low Reynolds number swimmers: An 
Oseen tensor approach, lattice Boltzmann and multipar- 
ticle collision dynamics. In Section IIII| each of these 
methods was used to model a very simple, linear swim- 
mer comprising three linked spheres. Analytic results are 
available for this system [3] and hence we were able to 
validate the approaches and identify the strengths and 
weaknesses of each method. For swimmers made up of 
a small number of spheres the Oseen tensor approach is 
very fast. However as the number of spheres increases, or 
as multiple swimmers are considered, lattice Boltzmann 
becomes more efficient. Moreover the Oseen tensor does 
not accurately take into account short range hydrody- 
namic interactions. Lattice Boltzmann is able to deal 
with spheres close to each other or to a surface, but at 
the expense of an increasingly fine lattice resolution. We 
found that multiparticle collision dynamics is in general a 
less useful approach as the intrinsic noise tends to domi- 
nate the results. This method will be of use when model- 
ing nanoscale systems where fluctuations are an intrinsic 
component of the physics. 

Subsequently, in Section IIV| we proposed a new three 
sphere swimmer, which has the advantage of being able 
to turn, and control its trajectory in a three dimensional 
manner. These ideas may be of use in the design and 
fabrication of artificial micros wimmers. Section [V] ex- 
tended the one dimensional three sphere swimmer aiming 
to search for the most efficient swimming strategy for a 
larger number of spheres. We found that the efficiency 
increases logarithmically with the sphere number. 

In Section IVT1 we looked at modeling a filament swim- 
mer that moves due to the propagation of waves along its 
length. Using multiparticle collisional dynamics we were 
able to reproduce the behavior observed experimentally. 

There are many directions in which it would be fruit- 
ful to pursue the simulations. We are currently consider- 
ing interactions between two or more swimmers, and it 
would be of interest to consider the effect of boundaries 
and obstacles on swimming behavior. Continuum hydro- 
dynamic theories have recently been pro posed to describe 
concentrated solutions of swimmers j34l |35[ . These have 
led to results very suggestive of swimming behavior but it 
is hard to pin down the phenomenological parameters in 
the equations of motion. Simulations of increasing num- 
bers of swimmers are needed to try to bridge the gap 
between the microscopic and continuum approaches. 

There is enormous scope presented by more realistic 
biological problems such as the chemotactic responses of 
bacteria or random tumbling during a swimming cycle. 
Moreover simulations of this type may be important in 
designing artificial microswimmers for applications as di- 
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verse as drug delivery or mixing fluids within microchan- nels. 
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